Non-Gaussian fluctuations in stochastic models with absorbing barriers 
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The dynamics of a one-dimensional stochastic model is studied in presence of an absorbing bound- 
ary. The distribution of fluctuations is analytically characterized within the generalized van Kampen 
expansion, accounting for higher order corrections beyond the conventional Gaussian approxima- 
tion. The theory is shown to successfully capture the non Gaussian traits of the sought distribution 
returning an excellent agreement with the simulations, for all times and arbitrarily close to the 
absorbing barrier. At large times, a compact analytical solution for the distribution of fluctuations 
is also obtained, bridging the gap with previous investigations, within the van Kampen picture and 
without resorting to alternative strategies, as elsewhere hypothesized. 
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Studying the dynamics of large ensemble of interacting 
entities is a fascinating field of investigations of broad ap- 
plied and fundamental interest. Atom and nuclei result 
for instance in extended systems of mutually interact- 
ing discrete elements. Similarly, proteins can be ideally 
imagined as a coherent sea of microscopic actors that 
e.g. densely populate the inside of the cells. Social and 
human communities are also examples of systems that 
display rather intricate dynamics, microscopically gov- 
erned by the very complex network of interlaced connec- 
tions among individuals. Surprisingly, and as follows a 
widespread observation of cross-disciplinary breath, reg- 
ular collective modes can spontaneously emerge at the 
macroscopic level, as resulting from the erratic dynamics 
of the microscopic discrete constituents [ll-Q- 

The time evolution of such inherent stochastic system 
is generally described in terms of a master equation 
a differential equation for the probability of observing 
the scrutinized system in a certain configuration at given 
time. In practically all cases of interest, solving the mas- 
ter equation proves a task of formidable complexity, and 
approximated strategies need to be implemented to gain 
analytical insight into the system being examined. The 
celebrated van Kampen expansion ^ represents a viable 
technique to enable one for analytical progress. It is cus- 
tomarily believed that such a method works efficiently 
well provided the system is defined in a open domain or, 
conversely, if it evolves sufficiently far from any existing 
boundaries As an emblematic example, when the 
system has to face an asymptotic extinction, thus evolv- 
ing towards an attractive absorbing state, the van Kam- 
pen perturbative scheme is only assumed appropriate for 
short times. Recently, and to eventually bypass these 



supposedly stringent limitations, a different perturbative 
approach has been pioneered in that suites for large 
times, when the system is feeling the absorbing barrier. 
Working within this generalized setting, and operating 
with reference to a paradigmatic model of systems with 
an absorbing state, the voter model, a closed expression 
for the distribution of fluctuations was obtained which 
agrees with direct simulations. The derivation rests how- 
ever on speculative grounds, which, despite the a poste- 
riori validation, seem to lack of a solid physical interpre- 
tation, fully justified from first principles. 

The purposes of this Letter are twofold. On the one 
side, and with reference to the same version of the voter 
model as considered in Q, we will analytically demon- 
strate that by extending the van Kampen expansion to 
include higher orders corrections (8|, beyond the classi- 
cal approximation, allows us to accurately reproduce the 
observed distribution of fluctuations at any time. Non 
Gaussian traits reflecting the presence of the absorbing 
barrier are nicely captured by the method, which proves 
therefore accurate also close to the boundary. Even more 
interesting, the van Kampen solution is shown to con- 
verge at late times to the distribution calculated in 0], 
this latter being hence explained within a sound and uni- 
versal descriptive picture. 

Let us start by introducing the stochastic discrete voter 
model. As in the spirit of [7| , we consider a system made 
of N elements in mutual interactions, possibly organized 
in different species. Label with Xi the elements of a 
specific species and with Xq all the other entities. The 
following chemical equations are proposed to rule the mi- 
croscopic dynamics: 
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The master equation which stems from the above sys- 
tem reads: 

d 



dt 



Pn{t) = (6--l)[r(n+l|n)P„(t)] 
+ {^+~\)[T{n~l\n)Pn{t)] 



(1) 



where P„ {t) is the probabihty of photographing the sys- 
tem at time i in a configuration with n individuals be- 
longing to the population of Xi and are the step op- 
erators . The transition rates are given by: 
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where the initial states are the right entries and the final 
states the left ones. As follows the above, n = is an 
absorbing state while n = N corresponds to a reflecting 
barrier. The van Kampen approach requires imposing: 
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where l/vN plays the role of a small parameter and 
paves the way to the perturbative expansion hereafter 
discussed. By inserting the working ansatz ([2]) into the 
master equation ([T]) , and hierarchically organizing the re- 
sulting terms with respect to their iV-dependence, one 
obtains at the first order the mean-field deterministic 
equation for the continuum concentration 0(t) (r being 
the rescaled time t/N), namely d4)/dT — —i'(j>, whose so- 
lution reads 0(r) = 4>q exp{—vT). Higher order contribu- 
tions results in a generalized Fokker-Planck equation for 

the new probability n(^, r) = P (^(/)(t) + ^/y/N, . By 
truncating the expansion at the second order yields the 
standard Fokker-Planck equation, which predicts Gaus- 
sian fluctuations. Allowing instead for higher order cor- 
rections, generates a cascade of terms whose relative 
weights are controlled by the finite size N. After a 
lengthy algebraic derivation one ends up with: 
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Formally, the positiveness of the probability n(-) is 
not guaranteed a priori under the generalized Fokker- 
Planck evolution, an observation that was made rigorous 
n,[l^. However, with reference to specific case studies 
it was shown that unphysical negative values are just 
occasionally attained by H(^, r), and punctually localized 
in the tails of the distribution. The phenomenon fades 
off when including a sufficiently large number of terms in 
the development. The adequacy of the prediction can be 
a posteriori evaluated via a direct comparison with the 
numerical experiments. 

To progress with the calculation, we set off to estimate 
the moments of the sought distribution n(^,r). Let us 
recall that the moment of order h is defined as 



Multiply both sides of the generalized Fokker-Planck 
equation by the factor and integrate over M in d^. 
A straightforward manipulation yields to: 
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where use has been made of the supposed regularity of 
the distribution H(^, r) to drop out the boundary terms 
resulting from integrating by parts 
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We therefore dispose of a closed system of first oder 
differential equations for the moments of the distribu- 
tion H(^,t). We can integrate it numerically and so 
estimate the quantities (^^), for all h, at any time r. 
The knowledge of the moments enables us to immedi- 
ately reconstruct the characteristic function, and so re- 
cover, upon Fourier transform inversion, the distribution 
H(^,t). The predicted profiles are displayed in Fig. [1] 
(solid line) for different times. A comparison is drawn 
with the outcome of direct stochastic simulations based 



on the exact Gillespie algorithm [12| (symbols), returning 
excellent agreement. The distribution of fluctuations dis- 
plays clear non-Gaussian traits. It gets in fact more and 
more skewed as time progresses, reflecting the non trivial 
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FIG. 1: The distribution of fluctuations at dis- 
tinct rescaled times r. The snapshots refer to: 
(a) r = 5, (b) r = 36, (c) r = 41 and (d) 
r = 50. The symbols stand for direct stochas- 
tic simulations. The solid lines represent the 
theoretical predictions as obtained within the 
generalized Fokker Planck scenario. We have 
in particular truncated the sums in the Fokker- 
Planck (O) to A: = 3 {l/N^^"^ corrections) and 
included 200 moments in the final estimates of 
the generating function. The dashed lines re- 
fer to the Gaussian solutions obtained work- 
ing within the van Kampen expansion at the, 
conventional, next to leading approximation 
terms). Here v = 0.01 and the dis- 
tributions are normalized so to have the maxi- 
mum equal to one. 



interplay with the absorbing boundary. Surprisingly, and 
at odds with what customarily believed, the van Kam- 
pen ansatz proves accurate well beyond the Gaussian ap- 
proximation that is often invoked to justify its intrinsic 
validity. As a corollary, it seems tempting to argue that 
the transformation ([2]) from discrete to continuum vari- 
ables is an exact one, and not just an approximation that 
presumably descends from the central limit theorem, as 
occasionally speculated. 

It is also very instructive to analyze the asymptotic fate 
of the distribution of fluctuations, as predicted within 
the realm of the van Kampen theory. Based on intuition, 
we expect that when time goes to infinity, the distribu- 
tion n(^,T) converges to a Dirac delta centered in zero. 
Indeed, plugging into the moments' equations ([3]), the 
asymptotic mean-field solution = 0, and looking for 
stationary solutions of the obtained system (i.e. setting 
the derivatives to zero), one readily gets {^^) — V/i, the 
moments of a delta function. However, for times large 
enough that ~ 0, but before the system has relaxed to 
its stationary state, the generalized Fokker-Planck equa- 
tion ^ reads: 
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where we have only retained the term in N dropping 
higher orders corrections. Perform now the scaling ^ — >■ 
C/ViN). The equation © can be cast in the form: 
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The large time distribution n(^',T) is therefore insen- 
sitive to the system size N and bears consequently uni- 
versal traits. Equation © can be solved analytically (see 
also [t'I) to give: 
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where /i(-) is the modified Bessel function of the first 
kind. For large r, recalling that I\{x) ~ xjl when x is 
small, one can approximate Eq. ([7]) as: 
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Operating with the rescaled variable which, it is 
worth emphasizing, emerges naturally within the van 
Kampem expansion, when the large time limit is being 
considered, it is equivalent to inserting into the govern- 
ing master equation the modified ansatz n = N <^ -\- 
This latter corresponds to the strategy adopted in Q for 
the specific choice a = 0. In other words, and interest- 
ingly enough, the expected fluctuations ^ are compara- 
ble to the discrete population size n, when the absorb- 
ing boundary is being approached. We have therefore 
recovered exactly the same solution as obtained in 
while working within the generalized, but conventional, 
van Kampen approach. The adequacy of ([8|) is challenged 
in Fig. [5] versus numerical simulations returning a per- 
fect quantitative agreement. Notice that different dis- 
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FIG. 2: The distribution of fluctuations at large times. Left panel: the distribution n(^, r) is plotted versus ^ at r = 388. 
Symbols refer to the simulations (A'^ = 1000), while the solid line stands for the (normalized) solution ^ after the change of 
variable ^' is performed. Right panel: the distribution 11 (^',r) is plotted as function of the rescaled Symbols 

refer to numerical simulations relative to distinct A*". In particular, — 500 (circles) and A'^ = 1000 (triangles). The solid line 
stands for the (normalized) solution ((8]). Here = 0.01 and the distributions are normalized to unit. 



tribution profiles recorded at distinct values N, nicely 
superpose when the rescaled fluctuations ^' is employed. 

As a side remark, we stress that the same conclusion 
can be drawn working in the Fourier space and oper- 
ating under analogous approximations. Retaining only 
1/y/N corrections in (j^l), assuming — > and perform- 
ing the scaling (ry'*) = {£,'^)N^^^, one immediately ob- 
tains the following N independent differential equations 
for the moments evolution: 

^(^") ^ -.MV') + fcll (2 - (9) 

Equations here obtained within the extended van 
Kampen scenario, could be also derived via the alter- 
native, supposedly distinct, approach discussed in Q. 
Equations (j9]) can be in fact straightforwardly deduced 
from Eq. ^ (see also Eq. (8) in [7|) following the same 
strategy for the evaluation of the moments as outlined 
above. 

In conclusion, the van Kampen approximation works 
effectively for all times, well beyond the Gaussian ap- 
proximation and in a regime where the presence of the 
absorbing boundary is definitely important. The method 
returns in fact the correct asymptotic solution (a delta 
function centered in the origin), but also converges to 
the large time solution calculated in f?\, which is there- 
fore contextualized within a general descriptive picture. 
At moderate times, after the Gaussian approximation 
has broken down, direct comparisons with the numer- 
ical experiments, as reported in this Letter, testify on 
the excellent predictive ability of the van Kampen the- 
ory. This is an important observation, that will certainly 
motivate using the van Kampen machinery beyond the 
limited domains of applications for which it was originally 
conceived, and/or later referred to Future investiga- 
tions will be targeted to extending the current analysis 



to cases where the absorbing boundary competes with a 
non trivial stable fixed point, as well as models for which 
the notion of space is explicitly accounted for . 
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